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The rheology and molecular structure of a model bitumen (Cooee bitu¬ 
men) under shear are investigated in the non-Newtonian regime using non¬ 
equilibrium molecular dynamics simulations. The shear viscosity, normal 
stress differences and pressure of the bitumen mixture are computed at dif¬ 
ferent shear rates and different temperatures. The model bitumen is shown 
to be a shear-thinning fluid at all temperatures. In addition, the Cooee 
model is able to reproduce experimental results showing the formation of 
nanoaggregates composed of stacks of flat aromatic molecules in bitumen. 
These nanoaggregates are immersed in a solvent of saturated hydrocarbon 
molecules. At a fixed temperature, the shear-shinning behavior is related to 
the inter- and intramolecular alignment of the solvent molecules, but also to 
the decrease of the average size of the nanoaggregates at high shear rates. 
The variation of the viscosity with temperature at different shear rates is 
also related to the size and relative composition of the nanoaggregates. The 
slight anisotropy of the whole sample due to the nanoaggregates is considered 
and quantified. Finally, the position of bitumen mixtures in the broad lit¬ 
erature of complex systems such as colloidal suspensions, polymer solutions 
and associating polymer networks is discussed. 
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I. INTRODUCTION 


Bitumen is one of the by-products of the refinery process of crude oil and a very 
viscous fluid. Its main application is as a binder in road pavement 1], because 
it combines many desirable properties: it is cheap, it adheres quite well to mineral 
filler particles and has suitable mechanical properties. Among bitumen’s mechanical 
properties, the shear viscosity is of special importance. To guarantee a cohesive 


pavement, the shear viscosity shou 
associated with a brittle behavior 


d not be too low, while high values are often 


2 ] and 3|] which facilitates the development of 
cracks and eventually leads to the deterioration of the pavement. Bitumen’s viscosity 
also has a great effect on the exploitation of bitumen reservoirs and transportation of 
bitumen. This is why the experimental literature reporting viscosity measurements 


of bitumen is so vast 


, 3—5]. Moreover, it is known that bitumen has a strong 


non-Newtonian behavior especially at low temperatures 


6-9]. 


The non-Newtonian 


behavior is characterized by a strain rate dependent shear viscosity and is believed 


to originate from the specific molecu 
particular asphaltene nanoaggregates 


ar structure of bitumen 10], containing in 

HI- 


Molecular dynamics (MD) is a good tool to study the molecular structure of 


bitumen and of the asphaltene nanoaggregates |12-|la], MD simulations can also 


quantify the mechanical properties of the material. Indeed, there are a few numerical 
studies relating the shear viscosity of bitumen in the limit of vanishing strain rate 
to its molecular and supramolecular structure 16dl9j. However, to our knowledge, 
there has been no MD study relating the non-Newtonian behavior of bitumen to its 
molecular structure, in particular to the nanoaggregate structure. That is the first 
aim of this article. A second aim is to compare the molecular structure of bitumen 
under shear to that of other complex systems such as colloidal suspensions 20j, 
polymer solutions 21|. and associative polymer networks { 22 ]. The hope is then 
to gain some insights on which structural characteristics of bitumen to focus on to 
model its rheology. 

In order to achieve these aims, we used non-equilibrium molecular dynamics sim¬ 
ulations with the SLLOD algorithm 23 and 124]. This algorithm is able to reproduce 


a shear Couette flow, even far from equilibrium 23] and provide the value of the 
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shear rate dependent viscosity. The SLLOD algorithm has been used to simulate 
shear flow of complex systems such as hyper branched polymers 25], associating 


polymer solutions 
bitumen 


26 ] and ionic liquids 22]. In this paper, we applied it to Cooee 


28], which is a four-component bitumen model describing a bitumen with 


generic mechanical properties 


17] and molecular structure 


18 


and 


29J. Different 


temperatures were investigated in order to extrapolate the results to ambient tem¬ 
perature. 

The paper is organised as follows. Section QT| contains information about the 
model bitumen used and the SLLOD algorithm. Results for the rheological prop¬ 
erties and the molecular structure of the sample are presented and discussed in 
Secs. HU and HVl Finally, Sec. M contains a summary and a discussion. 

II. SIMULATION DETAILS 


A. Cooee bitumen 


The chemical composition of bitumen is very complex 30]. However, the 
molecules in bitumen can be classified into dissolution classes. The SARA clas¬ 


sification 31] is one of the most common classifications and it distinguishes between 
Saturated hydrocarbons, resinous oil molecules, called Aromatic in the SARA 
scheme, Resin molecules and Asphaltene molecules. All of the last three molecule 
types are aromatic and the molecular weight of asphaltene molecules is on average 
higher than that of resin molecules, which is higher than the resinous oil molecu¬ 
lar weight. To obtain a simple model of a generic bitumen, one typical molecular 
structure for each class was chosen. The molecular structures chosen are shown 
in Fig. [D The reasons behind these specific choices are given in Ref. Q . The 
main system studied in this paper contains in weight fraction: 57.1% of saturated 
hydrocarbons, 7.6% of resinous oil molecules, 13.8% of resin molecules, and 21.5% of 
asphaltene molecules. The total number of united atom units in the system is 15570. 
The methyl (CH 3 ), methylene (CH 2 ), and methine (CH) groups are represented by 
the same united atom unit of molar mass 13.3 g-moD 1 and the sulfur atoms are 
represented by a united atom unit with a different molar mass 32 g-moD 1 . The 
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interaction potential between the united atom units contains four terms. The first 
term corresponds to an intermolecular potential. It is a Lennard-Jones poten¬ 
tial of the form Ulj(t ) = 4 e((<r/r ) 12 — (<r/r) 6 ) with parameters a = 3.75 A and 
e/k b = 75.4 K, where r is the interatomic distance and ks the Boltzmann constant. 
The method of shifted potential is used with a cutoff of 2.5 a. The three other terms 
of the interaction potential describe intramolecular interactions. They control the 
bond length between two connected particles, the angle between three consecutive 
particles, and the dihedral angle between four consecutive particles. They are given 
by the following expression: 


fdntra(r) 


\ Y k *( r v ~ 

bonds 

i Y kg (cos 9 

angles 



— cos 



5 

Y 5>cos"*. 

dihedrals 71=0 


( 1 ) 


The values of the parameters k s , lb , kg, 6q, and c n arc listed in a previous work m 
The simulations were performed at four different temperatures: 377 K, 452 K, 528 K, 
and 603 K and at constant density. For each temperature, the chosen density corre¬ 
sponds to an average equilibrium pressure equal to the atmospheric pressure. The 
four densities chosen are: 1007 kg-m -3 , 964 kg-m -3 , 894 kg-m -3 , and 833 kg-m -3 . 


B. SLLOD algorithm 


In this work, the molecular version of the SLLOD algorithm 23 and |25] is used 


to obtain a shear flow of the bitumen mixture. The equations of motion in that case 
are given by: 


Da 
P IOL 


D ■ Vu, 


Ha 


. rrin 


( 2 ) 

( 3 ) 


where r ia and p ia are the laboratory position and peculiar momentum of particle a 
in molecule i, respectively (peculiar means that the corresponding velocities are 
relative to the streaming/flow velocity), F, a is the force acting on particle a in 
molecule i, Vu is the velocity gradient tensor, chosen here such that the (x,y) 
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FIG. 1. (Color online) Structure of the four molecules in the ’’Cooee bitumen” model. 
Grey edges represent the carbon groups CH 3 , CH 2 and CH and yellow edges represent 
sulfur atoms. The “head” and “body” of the asphaltene molecule are shown. Numbers and 
arrows indicate bond-vectors used to quantify the nanoaggregate structure. Reproduced 
from Ref. (li|, (Copyright 2013 AIP). 


component is du x /dy = 7 and all other components are zero, r* = )C«=i m ia r ia/Mi 
is the position of the center of mass of molecule i, M t = , m,i a is the mass of 

molecule i, p* = Yla=i P*« * s the momentum of the center of mass of molecule i, Ni 
is the number of particles in molecule i, and ( is the Gaussian thermostat multiplier 
given by: 


c 


ES F, • Pi/Mi - 7 E :=1 PixPiy/Mi 


N„ 


ES P 1/Mi 


(4) 


where Fj = E^=i F JQ is the force acting on molecule % and N m the total number 
of molecules. The expression for the thermostat multiplier ( corresponds to the 
application of Gauss’ principle of least constraint, used to keep the molecular center 
of mass kinetic temperature constant. The algorithm used to propagate the SLLOD 
equations of motion is the operator splitting algorithm developed by Pan et al 32] 
and adapted for molecular SLLOD equations of motion. To our knowledge, it is 
the first time that the operator splitting algorithm has been adapted to molecu¬ 
lar systems. This implementation is described in appendix |A] The MD graphical 
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processing unit package RUMD 33] was used to perform the calculation in double 
precision. The time step is equal to 5t = 0.86 fs. Each simulation was equilibrated 
for a period of 1.7 ns and lasted thereafter 6.9 ns. Depending on the strain rate and 
the temperature, 8 or 16 independent initial configurations were used. The change 
in density from one initial configuration to the other is less than ICE 12 %. The 
averages and standard errors shown in this paper correspond to the averages and 
standard errors over these independent simulations. 


III. RHEOLOGY 

Four main rheological properties can be computed from the NEMD simulations 
under shear. They are the shear viscosity r), the first and second normal stress 
coefficients dq and d/ 2 , and the non-equilibrium molecular pressure P. They provide 
information about the Newtonian and non-Newtonian behavior of the fluid. 

These quantities are derived from the molecular stress tensor er, given by the 


Irving-Kirkwood formula 


34j: 


rr = 


1 

V 


Nm 

£ 

i =1 


Mi 


Nm N m 

y y ^, 

z=l j>i 


(5) 


where V is the volume of the system, (•) denotes a time average over the non¬ 
equilibrium steady state, r ij = r, — is the vector between the centers of mass of 
molecules i and j, F ;J - = F iajp is the intermolecular force acting on molecule i 

due to molecule j and F ia jp is the force acting on particle a in molecule i due to 
particle (3 in molecule j. 

A. Shear viscosity 


Due to the shear geometry enforced in the simulation, the shear viscosity is given 


by: 


V = 


g xv y ® 


xy i w yx 


27 


( 6 ) 


where 7 is the shear rate. The variation of the viscosity 7 with the shear rate 7 
is plotted in Fig. [2] for different temperatures. At all temperatures the viscosity is 
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decreasing with increasing shear rate and the Cooee bitumen is consequently a shear¬ 
thinning fluid. This is consistent with most experimental data on bitumen jfj, 8, and 
91. However, bitumen has a broad range of behavior with shear rate and can also 

n 

be shear-thickening |7J. Moreover, the higher the temperature, the higher the shear 
rate at which the fluid enters the Newtonian regime in the Cooee model. Thus, the 
Newtonian plateau is clearly visible at temperature 603 K and barely noticeable at 
temperature 377 K in the range of shear rates accessible to NEMD simulations. In 
other words, over the whole range of shear rates spanned, the change in viscosity is 


larger at low temperatures t 


lan at high temperatures. This is in agreement with 


the experimental literature js|, stating a more pronounced non-Newtonian behavior 
of bitumen at low temperatures. 

The zero shear rate viscosity p o can be determined using two empirical models 
to fit the data. The models are known as the Carreau-Yasuda model 35] and the 
Cross model 36]. The Carreau-Yasuda model suggests the following expression 
for the variation of viscosity with shear rate: p = Po/[l + (A7) 2 ] p , where A is a 
time constant and p a power law exponent. The Cross model suggests a different 
empirical expression for the variation of viscosity with shear rate: p = p^ + (r/ 0 — 
? ?oo)/[l + (K A f) m ], where p^ is the shear viscosity at infinite shear rate, K the 
consistency index, and m a power law exponent. The values of the parameters 
corresponding to the Carreau-Yasuda and Cross models are given in Tables [T] and HI] 
for all temperatures. The values of the zero shear rate viscosity p$ should be taken 
with care for the two lowest temperatures as the Newtonian regime is not clearly 
visible. 

The values of these different parameters will now be discussed. Both models 
obtain similar values for the viscosity at zero shear rate, which gives some confidence 
in the quality of the data. The experimental values for the zero shear rate viscosity 
at a given temperature depend a lot on the precise origin and age of the bitumen 


mixture studied. At around 377 K the experimental 
vary from 0.1 Pa-s [3] to more than 100 Pa-s (t] and 


values for the viscosity can 


16]. The value found in this 


work around 0.1 Pa-s at 377 K is consequently at the lowest end of the broad range 
obtained experimentally for different bitumen mixtures. 
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T [K] r/o [Pa-s] 

A [s] p 

377 

0.10 

2.9 x 10~ 8 0.29 

452 

0.022 

2.7 x 10” 8 0.21 

528 

0.0055 

2.6 x 10~ 8 0.13 

603 

0.0016 

6.7 x 10~ 9 0.09 


TABLE I. Values of parameters r/o, A, and p from the Carreau-Yasuda formula, for 
different temperatures T. 


The values of the time constant A in the Carreau-Yasuda model and of the con¬ 
sistency index K in the Cross model are also quite close to each other. They can 
be interpreted as the longest rotational relaxation time in the system j^]. The 
rotational relaxation time of the asphaltene molecules was estimated at tempera¬ 
ture 452 K and density 964 kg-m -3 at equilibrium in an earlier work jp|. It was 
evaluated as the time needed for the rotational correlation function to decay to 0.75 
and was found to be equal to r rot = 1.8 x 10~ 8 s. This value is very close to the 
value of the time constants A = 2.7 x 10“ 8 s and K = 3.4 x 10 -8 s found at the 
same temperature and density in this work from the Carreau-Yasuda and Cross fits, 
respectively. 

Finally, the Cross model does not seem to fit the data as well as the Carreau- 
Yasuda model for high shear rates, despite that the former model was designed 
to improve the fit at high shear rates by introducing a viscosity at infinite shear 
rate The fact that the viscosity keeps decreasing even at very high shear rates 
will be explained in terms of molecular aggregation in Sec. IIV Al 


To extrapolate the data to experimentally relevant temperatures around 300 K, 
we compare the variation of the zero shear rate viscosity with temperature to an Ar¬ 
rhenius behavior. The Arrhenius equation of the form r/ 0 (T) = r]°° exp(E/(kBT)), 
where rj°° is the viscosity at infinite temperature and E the activation energy, is 
usually valid for a range of temperatures well above the glass transition [Q. Ex¬ 
perimentally, the glass transition of bitumen depends on its exact composition and 
is typically around 250 K 38], which is well below the temperature range used in 
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FIG. 2. (Color online) Variation of shear viscosity 77 with the shear rate 7 for different 
temperatures. The viscosity values corresponding to the four lowest shear rates at tem¬ 
perature 603 K and 528 K were averaged over 16 independent simulations, whereas the 
other data points were averaged over 8 independent simulations. The value of the viscosity 
at the lowest shear rate and at temperature 603 K is not displayed because it is negative 
with an error bar larger than its absolute value. 


T [K] 7 ?0 [Pa-s] 

K[ s] 

7700 [Pa-s] 

m 

377 

0.13 

1.6 x 10~ 8 

0.0013 

0.79 

452 

0.040 

3.4 x 10~ 8 

0.00064 

0.55 

528 

0.0061 

4.2 x 10 -9 

0.00085 

0.70 

603 

0.0018 

6.2 x 10 _1 ° 

0.00040 

0.56 


TABLE II. Values of parameters 770 , K , 7^, and m from the Cross formula, for different 
temperatures T. 


this work. The Arrhenius equation has been shown to fit the variation of the exper¬ 


imental zero shear rate viscosity with tempera' 
a temperature range going from 273 to 363 K 


;ure quite well for some bitumen in 
8j. The variation with temperature 
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of the zero shear rate viscosity obtained with the Carreau-Yasuda and Cross fits 
is compared to an Arrhenius behavior in the inset of Fig. [3j Both sets of values 
follow an Arrhenius behavior. The parameters of the Arrhenius fit are found to be 
E/k B = 4073 K and r)°° = 2.27 x 10 -6 Pa-s for the Carreau-Yasuda set of data and 
E/k B = 4245 K and r)°° = 2.00 x 10 -6 Pa-s. Taking into account the data from 
the Carreau-Yasuda fit, the Arrhenius equation predicts a value for the viscosity 
of 2.5 Pa-s at 293 K. Again, it is at the lowest end of the broad range of experimen¬ 
tal values for the zero shear rate viscosity at this temperature. These values can 
range from 1 Pa-s 4] up to 10 4 — 10 6 Pa-s B 0 , y, 


and 


16]. 


The fact that the viscosity of the Cooee bitumen is quite low compared to exper¬ 
imental viscosities of bitumen at ambient pressure and temperature can be related 
to the observation that the density of the Cooee bitumen is itself quite different 
from experimental densities. At a temperature of 377 K, the density of the Cooee 
bitumen is 1007 kg-m -3 . It is to be compared with the thorough experiments carried 


out on Athabasca bitumen 


39] which find a linear dependency of bitumen density 


on temperature at ambient pressure. Extrapolating their results to 377 K leads to a 
bitumen density around 955 kg-m -3 , which is significantly lower than the density of 
the Cooee bitumen at the same pressure and temperature. This difference between 
the equation of state of our model and the equation of state obtained experimen¬ 
tally could be due to the specific type of asphaltene molecules that we chose. It 
has very short alkyl chains aside from the aromatic plane compared to the resin 


molecule. We know from a previous work 


18], that adding resin molecules while 


removing asphaltene molecules from the mixture has a tendency to lower the den¬ 
sity at constant pressure. Thus, considering asphaltene molecules with longer alkyl 
chains could lower the density of the model and could also change the viscosity. In 
the future, we plan to investigate this possible effect of the length of the alkyl chains 
on the density and viscosity of the mixture by comparing two mixtures with alkyl 
chains of different lengths. 

In an attempt to extrapolate not only the zero shear rate viscosity but also 
the shear rate dependent viscosity to ambient temperature, we tested the time- 
temperature superposition principle. This principle has been shown to be generally 
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valid exp erimentally 40] and is also discussed in another numerical model of bitu¬ 
men 41]. The reduced viscosity and reduced shear rate are defined as rj/a n and yax, 
respectively, where 

Vo(T, p) 


CLff — 


ax = cu 


^o(-^ref) Pref) 
-^refPref 


(7) 

( 8 ) 


v Tp ' 

The zero shear rate viscosity rj 0 is chosen from the Carreau-Yasuda fit, the reference 
temperature T ref and the reference density p ref are chosen to be 603 K and 833 kg-rn -3 
respectively, because the Newtonian regime at this state point is the most clearly 
defined. The procedure to obtain the reduced viscosity and reduced shear rate is 
described in Ref. 42] and its theoretical foundation is discussed in Ref. [ 35 1. In 
particular, the time temperature superposition principle can be shown to be valid if 
all the relaxation times in the system have the same temperature dependence 43]. 
Figure |3] shows the variation of the reduced viscosity rj/a^ with the reduced shear 
rate yax for the Cooee bitumen. The time temperature superposition principle is 
not so well satisfied for this model bitumen. 

To explain this fact, we tried first to take into account the colloidal nature of 
bitumen [ll] . Indeed, the Cooee bitumen mixture can be seen as a suspension 
of branched nanoaggregates composed of aromatic molecules and surrounded by a 


saturated hydrocarbon solvent 18. and [29]. In the experimental literature, bitumen 


mixtures are also sometimes treated as colloidal suspensions, because the typical 
size of a nanoaggregate is generally agreed to be around 2 nm 11 !]. hi the case of 
a dilute or semidilute suspension, the contribution of the suspended objects to the 
total viscosity can be well approximated by subtracting the pure solvent viscosity 
to the total viscosity. In the case of the Cooee bitumen studied here, the volume 
fraction of aromatic molecules involved in the nanoaggregates is 0 = 0.42, far from 
the dilute limit. We tried nevertheless to subtract the viscosity of the pure solvent at 
the same temperature and same density from the viscosity of the bitumen mixture, 
to see if the time temperature superposition principle was better satisfied. The 
reduced viscosity is now equal to (77 — rj s )/a' T) and the reduced shear rate to ya 0 , 
where rj s is the pure solvent viscosity at the same temperature and density and 
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FIG. 3. (Color online) Variation of the reduced viscosity r]/a ri with the reduced shear 
rate jax for different temperatures. This is a test of the time-temperature superposition 
principle. Inset: (Color online) Variation of the zero shear viscosity rjo, evaluated from 
both the Carreau-Yasuda and the Cross models, with the inverse temperature. The solid 
lines are Arrhenius fits to the data. 


where a' v and a 7 


are given by: 


a v = 


CLj 1 — CLjj 


Vo(T,p) -ri s (T,p) 

VO (TreU Pref) Vs ( -^ref* Pref) 

-^refPref 


Tp 


(9) 

( 10 ) 


The reduced viscosity is plotted versus the reduced shear rate in Fig. [U for the 
three highest temperatures. As expected, subtracting the solvent does not help 
the time superposition principle to be satisfied. Rather than a solvent effect, the 
failure of the time superposition principle will be related to the presence of at least 
two characteristic times related to the nanoaggregates and with a priori different 
temperature dependences in Sec. IIV Al The inset of Fig. [4] displays the variation of 
the viscosity of the pure solvent with shear rate for different temperatures. It can be 
seen that at the lowest temperature, the pure solvent never reaches the Newtonian 
plateau. It actually crystallizes. This crystallization is not seen in the bitumen 
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FIG. 4. (Color online) Variation of the reduced viscosity (77 — r] s )/a' 71 with the reduced 
shear rate 7 a' T at different temperatures. Inset: Variation of the viscosity r] s with the 


shear rate 7 at different temperatures for the pure solvent. 


mixture. There could be two main reasons behind this: the melting temperature of 
the mixture is lower than that of the pure solvent ; the bitumen mixture is blocked 
in a metastable liquid state or a glassy state during the time spans accessible to 
MD whereas the pure solvent is not. Regardless of the possible crystallization of 
the solvent at the lowest temperature, the fact that the bitumen mixture is so far 
from the dilute limit indicates strongly that aggregate-aggregate interactions are 
very important to understand the rheology of the mixture. 


B. Normal stress coefficients 

A fluid in the non-Newtonian regime is also characterized by shear-rate dependent 
normal stress coefficients. To first order in the shear rate, the two normal stress 
differences a xx — a yy and a yy — a zz are exactly zero because a shear strain does not 
cause any change in the diagonal components of the stress tensor. Different models 
and experiments on polymers show that the first departures of the normal stress 
differences from zero occur at second order in the strain rate [l(3]. Therefore, the 
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T[ K] 

a 

P 

377 

1.4 

1.4 

452 

1.2 

1.1 

528 

0.98 

0.81 

603 

0.96 

0.71 


TABLE III. Values of exponents a and /3 describing the power law decay of the normal 
stress coefficients Ti and \H 2 , respectively, for different temperatures T. 


normal stress coefficients 'Ll and T 2 are defined as: 'L 1 = (cr xx — cr yy )/ y 2 and d / 2 = 
(a yy — (Jzz)/i 2 - Thus, at low shear rates, the normal stress coefficients dy and T 2 
should plateau to a given value. At larger shear rates, the normal stress coefficients 
are known to decrease for polymers 10]. The normal stress coefficients dy and T 2 
obtained from the NEMD simulations of Cooee bitumen are shown in Fig. [5] (a) 
and (b), respectively. Only the results having a standard error lower than their 
absolute value are displayed. The first normal stress coefficient Ti decreases with 
increasing shear rate 7 , showing again that bitumen has reached a non-Newtonian 
regime in the range of shear rates spanned in the simulations. The second normal 
stress coefficient 'Id has been less studied, but is usually negative in experiments 
on polymer fluids [25 


and 


35]. The simulations presented in this paper also find 


a negative second normal stress coefficient T 2 . The absolute value of the second 
normal stress coefficient T 2 decreases with increasing shear rate, which is also a 
signature of the non-Newtonian regime. The decrease of the two normal stress 


25j: T, 


7 " 


coefficients in the non-Newtonian regime can be fitted to a power law 
and — T 2 ~ 7 -/3 for the coefficients dy and T 2 , respectively. The values of the 
exponents a and 6 are displayed in Table IIIII They are close to each other and 
decrease with increasing temperature, reaching a value lower than 1 at the highest 
temperatures. 

Finally, the normal stress ratio —\l/ 2 /\I / 1 is plotted versus shear rate in Fig. [5] (c). 

This ratio is known to be equal to 0.24 for linear polymer melts and to 0.3 for 

□ 

branched polymer melts pTj]. In the case of bitumen, it is higher, around 0.5, but 
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FIG. 5. (Color online) (a) and (b): Variation of the first \Fi and second ^2 normal stress 
coefficients with the shear rate 7 for different temperatures, respectively. In (b) — T 2 is 
plotted versus shear rate. The straight lines are power law fit to the data. Their exponents 
are given in Tabic HTT1 (c): Variation of the ratio —T 2 /T 1 with the shear rate 7 for the 
same temperatures. 


does not seem to depend on either shear rate or temperature. This value of the 


normal stress ratio —^ 2 /will be related to an anisotropy effect in Sec. IIV B II 


C. Pressure 


The pressure of the system is defined as minus a third of the trace of the molecular 
stress tensor given in Eq. (l5lh For simple fluids, the pressure P has been shown to 


have a power law dependency on shear rate 


44 and 


45j]: P — Pq + Irf, where P 0 is 


the pressure at equilibrium. For more complex fluids in particular for linear polymer 
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T[K] 

e 

377 

0.86 

452 

1.4 

528 

1.9 

603 

2.2 


TABLE IV. Values of the exponent e describing the power law increase of the normal 
isotropic pressure P for different temperatures T. 


melts, the pressure begins to c 


with a power law behavior 


46 


ecrease with increasing shear rate before increasing 
. The variation of the pressure with the shear rate 


is displayed in Fig. [6] for the Cooee bitumen model. The variation is found to be 
monotonic with the shear rate for the range of temperatures explored. It can be 
fitted quite well with a power law. The exponent of the power law is given for 
each temperature in Table IIV1 The observation that the pressure is monotonically 
increasing with the shear rate instead of having a clear minimum could be due to 
the fact that the molecules in bitumen are not all linear. Moreover, it is known that 
the pressure increases faster with shear rate for high generation dendrimers than for 
linear polymers jd6| . This fast increase can mask the local minimum of the pressure. 
The same mechanism could be at play in the bitumen mixture studied here. 

Figure |6] also shows that the density of the mixture at each temperature is chosen 
so that the equilibrium pressure P 0 is around the atmospheric pressure. 

Finally, for shear rates lower than 3 x 10 8 s -1 and for all temperatures, the 
pressure does not deviate significantly from its equilibrium value, enabling a direct 
comparison with experimental results, usually obtained at constant pressure. 


IV. MOLECULAR STRUCTURE 

As mentioned in Sec. Ill A! the Cooee bitumen model used in this paper contains 
four molecule types. In order of increasing molecular weight, these types are: do- 
cosane, resinous oil, resin and asphaltene. Of these, only the docosane molecules 
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FIG. 6 . (Color online) Variation of the normal isotropic pressure P with the shear rate 7 
for different temperatures. Inset: Variation of log(P — Pq ) with log( 7 ). Only data with 
error bars lower than the value of log(P — Pq) are shown. The dashed lines are straight line 
fits to this data. The solids lines in the main figure correspond to equation P = Pq + 67 e , 
where b and e were determined by the straight line fits to log(P — Po) versus log ( 7 ). The 
exponents e are given in Table HVl 


are not aromatic. This composition resembles the SARA classification 31] and is 
able to reproduce the characteristic supramolecular structure of bitumen in nanoag¬ 


gregates 


11]. The nanoaggregates are composed of the aromatic molecules aligned 


with respect to each other. An example of a nanoaggregate seen in the simula¬ 
tions is shown in Fig. [0 The docosane molecules can be seen as a solvent for the 
nanoaggregates. The aim of this section is to quantify the variation of the nanoag¬ 
gregate structure and of the inter- and intramolecular structures of each molecule 
type with shear rate and temperature as well as to compare it with the variation of 
the rheological properties obtained in Sec. IIIII 
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Asphaltene 


FIG. 7. (Color online) Snapshot of a linear nanoaggregate, obtained in equilibrium molec¬ 
ular dynamics. Asphaltene molecules are in blue, resin molecules in red and resinous oil 
molecules in green. Modified with permission from J. Chem. Phys. 141, 144308 (2014). 


A. Nanoaggregate size 


The definition of a nanoaggregate from a molecular point of view is described 


in detail in Refs. 


18 


and 


29] for equilibrium simulations. The nanoaggregates are 


composed of all the aromatic molecules and are branched. The branched structure 
arises from the asphaltene molecules having two flat parts, a head and a body, 
oriented in different directions. Each part can align to other flat aromatic molecules 
creating a branched structure jlij. In this work, we focus on the linear segments of 
these branched structures, because their size distribution has been studied in detail 
at equilibrium nr an earlier wor k Q, but tire presence of tire branched structure 
should be kept in mind. The definition of a linear segment sums up to the following 
rule: two aromatic molecules are nearest neighbors in the same nanoaggregate if they 
are “well-aligned” and “close enough”. Only the alignment of aromatic molecules 
to an asphaltene body and not to an asphaltene head is considered here, to extract 
purely linear structures. All molecules connected by this rule are part of the same 
linear nanoaggregate. The same rule was used in the NEMD simulations and the 
average number of aromatic molecules in a linear nanoaggregate could be quantified 
for different shear rates. The results are plotted in Fig. [8] (a) for all temperatures. 

The variation of the nanoaggregate size with the shear rate at a fixed temperature 
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(a) 


(b) 


FIG. 8 . (Color online) (a): Variation of the average number of aromatic molecules (IV mo i) 
in a linear nanoaggregate with the shear rate 7 for different temperatures. The open 
symbol and the dashed line represent the corresponding value at equilibrium and at tem¬ 
perature 452 K. (b): Probability P mo i to have a linear nanoaggregate containing N mo 1 
molecules versus N mo \ at temperature 452 K and for different shear rates from 6 x 10 7 
to 10 x 10 10 s 1 (the two lowest shear rates are omitted for the sake of clarity). The open 
symbols correspond to the same result at equilibrium. 


is first considered. At low shear rates, the size of the nanoaggregates is close to 
their size at equilibrium as is shown for temperature 452 K in Fig. [8] (a). The full 
distribution of the number of aromatic molecules in a linear nanoaggregate at low 
shear rates is also very close to the same distribution at equilibrium as can be seen in 
Fig. |8](b). At these low shear rates, the size_distribution has a biexponential shape, 


which we attributed, in a previous work 


291 ]. to longer nanoaggregates being less 


bent and containing more asphaltene molecules than smaller ones. After a plateau at 
low shear rates, the nanoaggregate size decreases quite sharply with increasing shear 
rate until it reaches a value close to 1, indicating that nearly all nanoaggregates are 
reduced to single molecules. In other words, the shear flow induces nanoaggregate 
rupture. Figure [8] (b) shows that at the highest shear rate and at temperature 452 K 
some longer nanoaggregates still remain but with a low probability. Interestingly, 
figure [8] (b) also shows that at high shear rates the size distribution is closer to 
a single exponential distribution. However, no noticeable variation of the relative 
composition of the nanoaggregates with their size was seen when the shear rate is 
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N m o, 

N mo , 

(a) 

(b) 



FIG. 9. (Color online) (a): Same as Fig. [3] (b) for a given shear rate of 6 x 10 7 s _1 and 
at different temperatures, (b): Variation of the fraction r\ of asphaltene molecules in a 
nanoaggregate with the number IV mo i of aromatic molecules in the same nanoaggregate, 
at different temperatures and at a fixed shear rate of 6 x 10' s . (c): Variation of the 
average number of aromatic molecules (A r mo i) in a linear nanoaggregate with the stress a xy 
for different temperatures. 


increased (not shown). 

The variation of the linear nanoaggregate size with temperature at a fixed shear 
rate is now considered. As can be seen in Fig. [3] (a), the nanoaggregate size at low 
shear rates decreases with increasing temperature. The size distribution of the linear 
nanoaggregates at a given low shear rate appears biexponential at all temperatures 
as can be seen in Fig. [9] (a). The cause of the decrease of the nanoaggregate size with 
temperature at low shear rates can be further investigated by looking at the rela¬ 
tive composition of the nanoaggregates for different temperatures. The fraction r a 
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of asphaltene molecule in a nanoaggregate is defined as the number of asphaltene 
molecules in this nanoaggregate divided by the total number of aromatic molecules 
in the same nanoaggregate. Figure [9] (b) shows the variation of the fraction t'a of 
asphaltene molecules with the nanoaggregate size for different temperatures and at 
a low shear rate of 6 x 10 7 s~h For linear nanoaggregates with a number of molecules 
larger than 3, the fraction ta of asphaltene molecules in the nanoaggregates increases 
with increasing temperature. Concurrently, the fraction of resin and resinous oil 
molecules in the nanoaggregates tends to decrease with increasing temperature in 
long nanoaggregates. It means that the overall decrease of the nanoaggregate size 
with temperature is mainly due to resin and resinous oil molecules being excluded 
from the aggregates at large temperatures. The asphaltene molecules stay aggre¬ 
gated even at large temperatures. It is probably an effect of the cohesive energy 
being higher in absolute value between two asphaltene molecules than between two 
resin or two resinous oil molecules. This is itself due the difference in the sizes 
of the aromatic structures. Asphaltene molecules have a larger aromatic structure 
than resin and resinous oil molecules, thus creating a larger 7T-stacking interaction 
between the aromatic planes of two asphaltene molecules than between two resin or 
resinous oil molecules. At high shear rates in Fig. [ 8 ] (a), it seems on the contrary that 
the nanoaggregate size increases with increasing temperature. It is due to the fact 
that the shear rate at which the nanoaggregates begin to break up is not the same 
at all temperatures. The higher the temperature, the lower the shear rate needed to 
break up the nanoaggregates. This apparent crossover between the curves giving the 
nanoaggregate size versus the shear rate at different temperatures can be removed 
if the average nanoaggregate size is plotted versus the shear stress a xy rather than 
the shear rate 7 . This can be seen in Fig. [9] (c). This figure also shows that above 
a certain shear stress, all curves collapse. Figure |9] (c) can be understood if one 
assumes that the size of the nanoaggregates depends on two factors: the thermal 
energy and the work done by the external shear force applied to the system. The 
thermal energy is directly proportional to temperature and the work done by the 
shear force can be estimated from the product of the shear stress a xy and the char¬ 
acteristic volume of an aromatic molecule in a nanoaggregate. At low shear stresses, 
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the work done by the shear force is negligible with respect to the thermal energy 
and the nanoaggregate size is determined by thermal energy alone. The nanoaggre¬ 
gates have then their equilibrium size. At high shear stresses, on the contrary, the 
thermal energy is negligible with respect to the work done by the shear force and 
the nanoaggregate size is determined by the shear stress a xy alone. The curves at 
different temperatures collapse at high shear stresses. At intermediate values of the 
shear stress, both factors play a role in the nanoaggregate size. In this intermediate 
regime, an increase in shear stress at constant temperature results in a monotonic 
decrease of the nanoaggregate size and an increase in temperature at constant shear 
stress results in a monotonic decrease of the nanoaggregate size. 

The connection between the rupture of the nanoaggregates around a given shear 
rate and the variation of the viscosity with shear rate and temperature quantified 
in Sec. m is now discussed. The overall decrease of the nanoaggregate size with 
increasing temperature at low shear rates contributes to the decrease of the viscosity 
observed with increasing temperature in Fig. [2j 


The rupture of the nanoaggregates at high shear rates is consistent with a decrease 
in the viscosity at these same shear rates instead of a plateauing around the viscosity 
at infinite shear rate 77 predicted by the Cross model. The Carreau-Yasuda model 
does not postulate the existence of a constant viscosity at infinite shear rate. It 
could be why it better fits the data than the Cross model, as was found in Sec. IIII1 


Moreover, the important change in the nanoaggregate size upon shearing could 
explain why the time-temperature superposition principle is not satisfied. Indeed, 
the time temperature superposition typically works when all relaxation times have 


the same temperature dependence 43]. In the case of the Cooee bitumen model, 


at least two main long relaxation times can be identified: the rotational relaxation 
time of the asphaltene molecules, close to the characteristic time A of the Carreau- 
Yasuda fit, and the inverse shear rate at which the nanoaggregates rupture. They 
have a priori different temperature dependences, which could explain why the time 
temperature superposition principle is not satisfied. 
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B. Molecular alignment 


The variation of the size of the nanoaggregates, a supramolecular structure spe¬ 
cific to bitumen, has been addressed in Sec. IIV Al This section (IIV HI) and the next 
(IIV C|) are devoted to more usual quantifications of the variation of the inter- and 
intramolecular structures with the shear rate and their link to the shear-thinning 
behavior reported in Sec. IIIII This section focuses on the alignment of the molecules 
in the flow direction, also known as form birefringence. 

Form birefringence is a departure from the isotropic equilibrium state of the fluid 
due to the applied velocity gradient. The anisotropy arises in this case from the 
alignment of the molecules in the flow direction 471 ]. The form birefringence is 
quantified using a molecular order tensor S m , defined for a given molecule type as: 


N„ 


s m = 


N„ 


E( u 


2=1 


u,-1 

3 


( 11 ) 


where N m is the total number of molecules of the considered type, u* is the unit 
vector denoting the principal direction of molecule i of the given type, I is the 
identity tensor, and (•) denotes a time average in the non-equilibrium steady state. 
The principal direction Uj of molecule i is the unit eigenvector corresponding to the 


largest eigenvalue of the gyration tensor R^, defined as: 


p2 _ 1 

9 2 Ni 


Ni Ni 


- E E< r « - r /5) ® ( r » - r «)> 


( 12 ) 


a =1 /3=1 


where N t is the number of atoms in molecule i and r Q is the position of atom a in 
molecule i. Once the molecular order tensor S m is determined for each molecule type, 
two quantities derived from it are of interest: the molecular order parameter S m and 
the molecular alignment angle Xm- The molecular order parameter S m is defined 
as 3/2 times the largest eigenvalue of the molecular order tensor. It is equal to 0 
when the molecules are randomly oriented and to 1 when all molecules are perfectly 
aligned. The molecular alignment angle Xm is equal to arctan(e y /e a; ), where e = 
(e x , e y , e z ) is the eigenvector corresponding to the largest eigenvalue of the molecular 
order tensor S m . The molecular alignment angle Xm represents the angle between 
the average molecular orientation e and the shear direction (^-direction in this case) 
in the xy- plane. 
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1. Molecular order parameter 


The variation of the molecular order parameter S m with the shear rate 7 is dis¬ 
played in Fig. [10] (a) for a temperature of 452 K and for the different molecule 
types. The curves are quite different from one molecule type to another. For do- 
cosane molecules, the molecular order parameter S m tends to zero in the limit of 
vanishing shear rate. This behavior is expected, since the docosane molecules are 
isotropically distributed at equilibrium. The molecular order parameter for docosane 
increases monotonically to reach the value 0.55 at the highest shear rate simulated, 
showing an alignment of docosane molecules in the shear direction. This is con¬ 


sistent with the shear thinning behavior globally observed for Cooee bitumen 


48]. 


For asphaltene, resin and resinous oil molecules the curve is qualitatively different. 
At low shear rates, until a critical shear rate of approximately 7 = 2 x 10 9 s -1 at 
temperature 452 K is reached, the molecular order parameter of the three aromatic 
types is roughly constant and non-zero. Above this critical shear rate, the molecu¬ 
lar order parameter begins to increase with the shear rate. The small finite value 
reached by the molecular order parameter at low shear rates is equal to the one 
obtained from equilibrium molecular dynamics simulations, as can also be seen in 
Fig. [10] This value is higher for asphaltene molecules than for resinous oil molecules 
and for resin molecules. This unusual behavior of the molecular order parameter is 
easily explained if the presence of nanoaggregates at low shear rates is considered. 
Indeed, the long nanoaggregates tend to align with respect to each other, due to 
steric hindrance, creating a slight anisotropy in the system. This anisotropy is also 
present at equilibrium. At equilibrium, it was previously determined [29|] that long 
nanoaggregates contain more asphaltene molecules than resinous oil molecules and 
that the resin content is the smallest. This explains why the anisotropy is larger 
for asphaltene molecules than for resinous oil molecules and the smallest for resin 
molecules. According to Fig. [ 8 ] (a), at temperature 452 K, the nanoaggregates begin 
to break up at a shear rate very close to the critical shear rate j c = 2 x 10 9 s -1 . This 
indicates that the increase of the molecular order parameter for shear rates higher 
than 7 C is due to the alignment of the single molecules in the flow direction. Again, 
this is consistent with the shear thinning behavior observed at high shear rates. 
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FIG. 10. (Color online) (a): Variation of the molecular order parameter S m with shear 
rate 7 for different molecule types and for a temperature of 452 K. The open symbols 
represent the corresponding values at equilibrium, (b): Variation of the ratio r = s ro t/s C onf 
comparing the standard deviation due to the anisotropy of the sample to the standard 
deviation due to different samples with the shear rate 7 at a temperature of 452 K. The 
error bars are calculated for each type of standard deviations as the absolute value of the 
difference between the standard deviation corresponding to half of the data set and the 
standard deviation corresponding to the other half of the data set. The relative error on 
the ratio r is then calculated as the sum of the relative errors on s rot and s con f. The dashed 
line has equation r = 1 and corresponds to the anisotropy of the sample being the main 
source of variation between one initial configuration and another. 


Moreover, the existence of this intrinsic anisotropy in the sample could explain why 
the value of the normal stress ratio — is different from what is found for linear 

or branched polymer melts. Indeed, the normal stress ratio — is a measure 

of the anisotropy between the xx- and zz- directions on the one hand and the xx- 
and yy- directions on the other hand, which can be influenced by an ordering of the 
system. 

With the MD simulations which we carried out, we cannot know whether the 
global anisotropy of the system is a finite size effect or if it is also seen in macro¬ 


scopic samples. In this last case, bitumen could resemble the nematic p 


rase of a 


discotic liquid crystal, formed by the stacking of disk-shaped molecules 49]. Ex¬ 


perimentally, it is known that there is a change of the supramolecular structure of 
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bitumen with the concentration of asphaltene molecules. At low concentrations, 
small and isotropically distributed nanoaggregates can be seen. At larger concen¬ 
trations, the nanoaggregates gather in clusters. At even larger concentrations, the 
clusters are all connected to form a network [ll]. This last stage could be repre¬ 
sentative of the simulations carried out in this work. The global anisotropy of the 
system was evaluated using a molecular order tensor Q defined in a very similar way 
than the molecular order tensor S m , specified for each molecule type in Eq. (fTTTh 
The only differences are that the average is done over all aromatic molecules instead 
of over one particular type and that each molecule orientation is taken to be the unit 
vector normal to the molecule plane instead of the principal direction determined by 
the tensor of gyration. The molecular order parameter Q is then 3/2 of the largest 
eigenvalue of the molecular order tensor Q. It is found to be equal to Q = 0.12±0.01 
at equilibrium 29]. This value of Q is sufficiently small, so that the shear viscosity 
can be treated as as an isotropic scalar quantity. 

To further investigate the effect of the anisotropy of the sample on the viscosity, 
simulations were performed in the same shearing geometry as before and with the 
same starting configurations except that the y- and z- coordinates are exchanged. 
The simulations were carried out at a temperature of 452 K for the same shear 
rates and the same density. The shear viscosity was analyzed and the following 
comparison was made. Let ry and 77 ' denote the viscosities of the system starting 
with configuration i and the corresponding rotated configuration, respectively. Let 
also fii be the average ( 77 ; + r]' i )/2. The standard deviation s rot due to the anisotropy 
of the configuration is quantified as: 


■Tot 


\ 


N, 


conf 


-^conf 

M “ 7i) 2 . 


(13) 


i= 1 


where N con f = 8 is the total number of initial configurations considered. A similar 
standard deviation s con f comparing different initial configurations and not an initial 
configuration with a rotation of itself is defined as: 


^conf 




(14) 


where i and j are two random different configurations, chosen from the 8 different 
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FIG. 11. (Color online) Variation of the molecular order parameter S m with shear rate 7 


at different temperatures for asphaltene (a), resin (b), resinous oil (c) and docosane (d) 


molecules. 


configurations available. Exactly iV con f pairs (i, j) are considered for the sake of 
symmetry with the definition of s rot . The variation of the ratio r = s r ot/sconf with 
the shear rate is displayed in Fig. [TO] (b). The ratio of the two standard deviations 
is around 1 for all shear rates indicating that the main source of difference between 
two configurations is the anisotropy of the sample. The values of the viscosity 
given in Fig. [2] should therefore be understood as an average viscosity over different 
orientations. 

Figure [TT] shows the variation of the molecular order parameter of the four dif¬ 
ferent molecule types with shear rate for different temperatures. From Fig. [IT] it 
is quite clear that for all molecule types, the molecular order parameter at a given 
high shear rate increases with decreasing temperature. This is due to the fact that 





















for a given shear rate at low temperatures the stress felt by the molecules is higher 
and the alignment in the direction of the flow is more pronounced than at high 
temperatures. At low shear rates, the behavior is quite different from one molecule 
type to another. For docosane molecules, the molecular order parameter S m (D) 
decreases with increasing temperature, for the same reason as at high shear rates. 
For resin and resinous oil molecules the molecular order parameter also decreases 
with increasing temperature at low shear rates. This is probably because as the 
temperature increases, the resin and resinous oil molecules are more often present 
as single molecules than in long nanoaggregates (see Fig [8] (d)), which reduces their 
global ordering. In contrast, asphaltene molecules stay in long nanoaggregates even 
at high temperatures and their molecular order parameter is nearly independent of 
temperature at low shear rates. 

An interesting feature can be noticed for the molecular order parameter S m (A) of 
asphaltene molecules at high temperatures. Looking closely at the data of Fig, fill (a) 
at temperature 603 K, one can note three different regimes with increasing shear 
rate: at low shear rates S m (A) is mostly constant, at intermediate shear rates S m (A) 
slightly decreases with the shear rate, and finally at large shear rates it increases. 
The intermediate shear rate range corresponds roughly to [3 x 10 s , 6 x 10 9 ] s -1 
for temperature 603 K. At these shear rates, the average number of molecules in 
the nanoaggregates is still very close to its equilibrium value as can be seen in 
Fig. [8] (a). Consequently, the decrease in the molecular order parameter of the 
asphaltene molecules in this intermediate shear rate range can be interpreted in 
the following way: the nanoaggregates are still well formed, but the shear flow 
is strong enough to disturb the alignment of the nanoaggregates with respect to 
each other causing the order parameter to locally drop. At higher shear rates the 
nanoaggregates break, the molecules align in the shear direction and the molecular 
order parameter increases again. This intermediate drop of the molecular order 
parameter S m (A) of asphaltene molecules is also noticeable for temperature 528 K, 
but not for lower temperatures. It might be because the nanoaggregates break before 
the shear flow is strong enough to modify the alignment of the nanoaggregates 
with respect to each other. The direct consequence of this reorganization of the 
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nanoaggregates with respect to each other on the behavior of the viscosity with the 
shear rate is unclear. In total, the fluid is shear-thinning in this range of shear rates, 
as can be seen in Fig. [2] This global behavior could be mainly due to the alignment 
of the docosane molecules in the shear direction, the order parameter of the resin 
and resinous oil molecules being nearly constant in this range. 


2. Molecular alignment angle 


The variation of the molecular alignment angle Xm is plotted for each molecule 
type versus the shear rate 7 at temperature 452 K in Fig. [12] (a). The typical behav¬ 
ior for a polymer solution under shear is the following: in the Newtonian regime, the 
molecular alignment angle is around 45 degrees and in the non-Newtonian regime, 
the molecular alignment angle decreases until reaching zero. The value of 45 de¬ 
grees is expected in the Newtonian regime, due to the nonuniform distribution of 


angular velocity of the molecules in the shear flow 


25]. At larger shear rates, the 


molecules align to the shear flow resulting in a vanishing molecular alignment angle. 
This general behavior is roughly seen for all molecule types in Fig. [12] (a), though 
the error bars are very large at low shear rates. The large error bars at low shear 
rates can be explained by the fact that at these shear rates the nanoaggregates are 
well formed and tend to align with respect to each other, giving rise to the slight 
anisotropy of the whole sample. The alignment of the nanoaggregates with respect 
to each other causes the set of the angles with the shear direction to be ill-spanned 
compared to an isotropic material. The error bars decrease when the nanoaggregates 
are ruptured. This happens at quite large shear rates, when the single molecules 
begin to align in the flow direction, so that the angle of 45 degrees is not clearly 
visible. Furthermore, at high shear rates, the values of the molecular alignment 
angle y rn for different molecule types are ordered according to the molecular weight, 
the molecules with the highest weight having the lowest alignment angle. It is ex¬ 
pected as the difference in velocity from one end of the molecule to the other is 
larger for large molecules than for small molecules, giving rise to an alignment in 
the flow direction at lower shear rates. 
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Finally, Figs. [121(b) and (c) show the molecular alignment angle of asphaltene and 
docosane molecules, respectively, at different temperatures. The molecular align¬ 
ment angle of resin and resinous oil molecules at different temperatures is not dis¬ 
played for the sake of readability and resembles that of the asphaltene molecules at 
different temperatures. For both asphaltene and docosane molecules, the lower the 
temperature, the closer to zero the angle to the flow direction for a given high shear 
rate. This is consistent with the molecular order parameter being higher for low 
temperatures than for high temperatures at a given high shear rate. The error bars 
on the alignment angle of asphaltene molecules are large up to very high shear rates 
for high temperatures, consistent with the nanoaggregates being maintained up to 
these high shear rates. Due to the large error bars, the existence of the intermediate 
range of shear rates identified in Fig. [TO] (a) and corresponding to the shear flow 
disturbing the alignment of the nanoaggregates with respect to each other, does not 
have a clear effect on the alignment angle of the asphaltene molecules in the flow 
direction. 


C. Intramolecular structure 

The intramolecular structure can be characterized by the radius of gyration of 
each molecule type. The radius of gyration is defined here as (Rg) = trace(R^), 
where R~ is the tensor of gyration given in Eq. (1121) and (•) is an average over 
molecules of the same type and over time. Figure [13] shows the variation of the radius 
of gyration of each molecule type with the shear rate for different temperatures. 
The variation of the radius of gyration with temperature and with shear rate is 
very different from one molecule type to another. For asphaltene and resinous oil 
molecules, Figs. H3] (a) and (c) respectively, the radius of gyration barely changes 
with the shear rate. It is expected as nearly all carbons in these molecules are 
aromatic creating a stiff intramolecular structure. However, at the largest shear 
rates the radius of gyration increases slightly indicating a stretching of the bonds. 
For both molecule types the radius of gyration increases with increasing temperature, 
indicating a stretch of the stiff bonds due to a higher velocity of each bead. 
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FIG. 12. (Color online) Variation of the molecular alignment angle Xm with shear rate 7 


for three cases: at temperature 452 K and for different molecule types (a), for asphaltene 


molecules at different temperatures (b), and for docosane molecules at different tempera¬ 


tures (c). 


For docosane and resin molecules, Figs. [13] (b) and (d) respectively, the trends 
are very different. Both molecule types can possess a maximum of their radius of 
gyration at a given shear rate. For linear molecules such as docosane molecules, this 
behavior is known 


46| 


and 


51] and is attributed to a stretching of the bonds in the 


flow direction for intermediate shear rates and a coiling up of the whole molecule at 
very high shear rates when each molecule tumbles more often. The maximum in the 
radius of gyration of linear molecules is usually associated with a clear minimum of 


the pressure of the system 
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and 


5,1]. It is not seen in the case of the Cooee bitumen 


model, probably because other molecules play a part in the value of the pressure 
as mentioned in Sec. IIII Cl We attribute the maximum seen for resin molecules at 
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FIG. 13. (Color online) Variation of the radius of gyration with the shear rate for asphal¬ 


tene (a), resin (b), resinous oil (c) and docosane (d) molecules at different temperatures. 


a given shear rate to a similar effect, because resin molecules have some long linear 
alkyl chains aside from the aromatic plane, in contrast to resinous oil and asphaltene 
molecules. Finally, for both docosane and resin molecules, the radius of gyration 
decreases with increasing temperature. This could also be an effect of the larger 
velocity of each bead, causing coiling of the molecules when the bonds are less stiff 
than in pure aromatic molecules. 

Another informative quantity which can be defined to describe the intramolecular 
structure is the bond order tensor Shonci- It is defined for docosane molecules as: 


N b -1 


S bond - N T 


a =1 


-i 

3 


(15) 


where N .& is the number of bonds in a docosane molecule, (•) is an average over 
docosane molecules and over time, and v a is the unit vector between neighboring 
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atoms a and a + 1 given as: 


U+i ~ r a 

v« = i- 

|r Q +i - r Q 

In a similar way as for the molecular order tensor S m , a bond order parameter S'bond 
and a bond alignment angle y b ond can be defined from the bond order tensor S b 0 nd- 
These two quantities measure the extent of the alignment of intramolecular bonds 
to the shear direction and the angle between them. Figures [M] (a) and (b) show the 
variation of the bond order parameter S’bond and the bond alignment angle Xbond of 
docosane molecules with the shear rate for different temperatures. Both quantities 
have similar variations with shear rate and temperature as their molecular counter¬ 
parts. As the shear rate increases, the bonds inside the docosane molecules become 
more and more aligned in the direction of the shear. This effect gets stronger as 
the temperature decreases. It is also consistent with the shear-thinning behavior 
observed in Sec. IIIIA1 

The bond order tensor associated with aromatic molecules is less well defined, 
because the bonds form rings and are not clearly oriented from one end of the 
molecule to the other. That is why it is not calculated here. However, as shown by 
the variation of the radius of gyration of the aromatic molecules with the shear rate, 
the stretching of the molecules is less pronounced for stiff aromatic molecules than 
for linear ones. Changes in the intramolecular structure of the aromatic molecules 
probably contribute very little to the shear-thinning behavior generally observed. 



V. SUMMARY AND DISCUSSION 


Non-equilibrium molecular dynamics simulations under shear have been per¬ 
formed for the first time on a bitumen mixture. They reveal both the non-Newtonian 
regime of such mixtures and their microscopic structure. The non-Newtonian regime 
is characterized by a shear rate dependent viscosity, decreasing normal stress coef¬ 
ficients, and a monotonically increasing pressure. The Cooee bitumen is found to 
be shear-thinning at all temperatures investigated in line with most experimental 
results on bitumen jg]. The shear-thinning behavior is in agreement with the varia¬ 
tion of some characteristics of the inter- and intramolecular structure with the shear 
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(a) (b) 


FIG. 14. (Color online) Variation of the bond order parameter Sbond (a) and the bond 
alignment angle Xbond (b) with shear rate 7 for docosane molecules at different tempera¬ 


tures. 


rate: 


• The saturated hydrocarbon molecules, which can be seen as solvent to the 
nanoaggregates, align with respect to each other and align their bonds in the 
shear direction, at all shear rates. 

• At high shear rates, the nanoaggregates disintegrate and the single aromatic 
molecules align in the shear direction. 

Some characteristics of the inter- and intramolecular structure of the Cooee bitumen 
do not change with the shear rate and do not seem to prevent the shear-thinning 
behavior from being established. They are: 

• Until a critical shear rate, dependent on the temperature, the nanoaggregate 
size distribution is unchanged with increasing shear rate. 

• The intramolecular structure of the aromatic molecules, composing the nanoag¬ 
gregates, is mainly constant up to quite high shear rates. 

Finally at high temperatures, there is a domain of intermediate shear rates, where 
the nanoaggregates still have their equilibrium size, but where the molecular order 
parameter of asphaltene molecules decreases with increasing shear rate. The de¬ 
crease of the molecular order parameter indicates a perturbation of the alignment 
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of the nanoaggregates with respect to each other, due to steric hindrance, by the 
shear flow. The effect of this local reorganization on the variation of the viscosity 
with the shear rate is unclear and could be masked by the inter- and intramolecular 
alignment of the solvent molecules. 

In addition, the variation of the rheological properties and of the corresponding 
molecular structure with temperature was studied. As expected, the viscosity de¬ 
creases with increasing temperature, following an Arrhenius law. The nanoaggregate 
size decreases with increasing temperature, which also contributes to the decrease 
of viscosity with temperature. The increase in viscosity at low temperatures and at 
a given shear rate explains some of the characteristics of the molecular structure: 

• The solvent molecules begin to align in the shear direction at lower shear rates 
at low temperatures than at high temperatures. 

• The nanoaggregates break up at a lower shear rate at low temperatures than 
at high temperatures. This is also due to the nanoaggregates being longer at 
low temperatures than at high temperatures. 


The study of the structure and rheology of the bitumen mixtures at different 
temperatures enables us to compare the behavior of a bitumen mixture with other 
complex systems such as colloidal suspensions, polymer solutions and associative 
polymer networks. The comparison with these complex systems is relevant because 
they have similarities with a bitumen mixture. Indeed, the Cooee bitumen mixture 
can be seen as a solution of patt.y flexible —gates of various s.ses Q, with 
a branched structure [18j, in a saturated hydrocarbon solvent, as was already men¬ 
tioned in Sec. IIV A[ Furthermore, the nanoaggregate size depends on temperature 
and shear rate as was seen in Sec. m 

It is known that both the linear rheological properties such as the zero shear rate 
viscosity and the non-linear rheological properties of colloidal suspensions and poly¬ 
mer solutions depend on the size, the shape and the polydispersity of the dispersed 


objects 


0,E, 


and 


53]. In the limit of a dilute solution of such objects, a lot is known 


on the relation between the intrinsic viscosity [ 77 ], defined by 770 = r] 0s + [ 77 ] 0 + o( 0 ), 
where 7]q s is the zero shear rate viscosity of the solvent and 0 the volume fraction 
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of objects, and the shape of the colloids or the topology of the polymers [53]. The 
shear-thinning behavior at high shear rates can also be related to the topology of the 
polymers, though for complicated topology such as dendrimers, quantitative models 
are not available yet to predict the exponent of the decay of the viscosity in the 


limit of infinite shear rate 


21 J. At high volume fractions of objects, a higher order 


expansion of these linear and non-linear rheological properties with the volume frac¬ 
tion is needed, because the interactions between the colloids or polymers cannot be 


neglected 


20 ], Bitumen mixtures typically fall in this range, which makes a direct 


and quantitative comparison with what is known of the rheological properties of 
dilute and semi-dilute polymer solutions of different topologies difficult. 


Another characteristic of bitumen which makes a direct comparison with col¬ 
loidal suspensions or polymer solutions difficult is the fact that nanoaggregates are 
maintained by weak interactions and their size depends on temperature and shear 


rate. T 


cal 


lis is usually not the case for colloids and for polymer solutions. Theoreti- 


m 


55 


57] and numerical 


26 


58 


62 ] studies have been carried out on the linear 


and non-linear rheological properties of solutions of associative polymers, where the 
size of the assemblies also depends on temperature and shear rate. There is for 


example 


571 . and 


a growing literature on the rheo 


60] and of wormlike micelles 


58 


ogv 


of telechelic polymer networks 22 


59, 61, 


55- 


and 


62] , However, the knowledge 


in this field is still scarce and very specific to the precise topology of self-assembly 


studied 


261 ]. It is consequently not so easily transferable to bitumen. 


Nevertheless, the current knowledge on the rheology of colloidal suspensions, 
polymer solutions and associative polymer solutions defines a framework which can 
be used to better understand bitumen. The size, the polydispersity and the topology 
of the nanoaggregates as well as the interaction between the nanoaggregates are 
expected to be the important parameters controlling the rheological properties of 
bitumen. We hope that the numerical study carried out in this paper can be used 
as a step towards an analytical modeling of the rheological properties of bitumen 
based on these parameters. 
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Appendix A: Operating splitting algorithm for molecular SLLOD 


We have implemented a molecular version of the operator-splitting algorithm 


described in Ref 


32 


To our knowledge, it is the first time that this algorithm is 


adapted to molecular systems. The implementation involves a two-step process: 
(1) the equations of motion of the molecular centers of mass are solved; (2) the 
(relatively simple) equations of motion for the particles’ positions and velocities 
relative to the molecular centers of mass are solved and combined with the solutions 
for centers of mass in order to update the atomic coordinates and velocities. Note 
that here velocities are the so-called peculiar velocities, i.e. velocities relative to the 
streaming velocity. 

The operator splitting method is an approximation based on a factorization of 
time propagator U generated by the system’s Liouville operator iL. The factor¬ 
ization is denoted AB 1 B 2 B 1 A where A is the operator associated with evolution 
of the coordinates at fixed momenta, B x is that associated with evolution of the 
momenta due to the flow only, and B 2 is that associated with the evolution of mo¬ 
menta in the absence of flow, at fixed particle coordinates. Exact solutions of each 
set of equations of motion are given by Pan et Note that A and Bi appear 
twice and thus should be integrated over a half time-step each time, although the 
first instance of A can be “wrapped” around and put at the end, at the expense 
of having the positions and velocities and positions a half-step out of sync. The 
isokinetic thermostat is included in each factor, so that the kinetic energy is exactly 
conserved (up to round-off error). 


38 



For the molecular case the same factorization is made and the above two step 
procedure is applied for each factor in the approximate propagator. The equations 
of motion for operator A are, for our choice of simple shear deformation (where i is 
a molecular index, a is a particle index within a molecule, and quantities with only 
the former are center-of-mass quantities): 

= — + 7 UiX = Via + 7 Vi^, (Al) 

Wlia. 

P ia = 0 (A2) 


For later convenience we introduce the atomic and molecular velocities Vj a and v,; 
respectively. Focusing on the positions, we multiply by the mass fraction of 

the particle a in molecule i and sum over a to get the equation of motion for the 
center of mass position: 


Pi . . . 

K = — + IViX = Vi + WiX. 


(A3) 


This is identical to the position part of the A operator in Ref. [32j (re-interpreting 
the index i in their equations as a molecular index), and has the same solution for a 
time increment At (although incorrectly given in that paper; the correct expressions 


can be found in Ref. 


631): 


Xi[t + At) = Xi(t) + A t(p xi (t)/Mi + 7i/i(t)) + (l/2Mi)At 2 ^p yi (t) (A4) 

Vi(t + At) = Ui(t) + A t(pyi(t)/Mi) (A5) 

Zi(t + At) = Zi(t) + A t(p zi (t)/Mi) (A6) 


Knowing the evolution of the center of mass position during the time step At, we 
can now consider the equations of motion for the relative coordinates of the particles 
within the molecule by subtracting Eq. (I A 31) from Eq. (1 A If) . giving 


d 

dt 


(r ia 




P ia/m ia - Pi /Mi = V ia - Vi, 


(AT) 


where the right side is the particle’s velocity relative to the molecule’s center of 
mass and is constant for the A equations of motion. The solution is therefore (using 
velocities instead of momenta) 


Yiait + At) - Ti(t + At) = r ia {t) ~ Yi{t) + A t(v ia (t) - Vi (f)), (A8) 
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which, upon substituting Eqs. (1A4j) - IA6l gives 


x ia (t + At) = x ia (t) + A t(v xia (t) + 7 Vi(t)) + l/2At 2 'yv yi (t) (A9) 

Viait + At) = y ia (t) + A t(vyia(t)) (A10) 

Z ia (t + At) = Zia(t) + A t(v zia (t)) (All) 

The same procedure is applied to the equations of motion for the Bj operator, 

which for the atomic momenta (at constant positions) are: 

dpi, 


mi 


mi 


dt =- 7 MT^ -0il M 


Pi 


(A12) 


Here aq is the part of the Gaussian thermostat multiplier Q associated with conserv¬ 
ing the kinetic energy during this part of the integration; it has the expression 

Ei IPxiPyi/Mi 


ot.\ — 


(A 13) 


E,pf/v, 

Note that the thermostat term involves only molecular quantities, that is, it is a 
so-called molecular thermostat which conserves the sum of the translational kinetic 
energies of the molecules ESumming Eq. (1A12I) over a gives 


dpi 

dt 


which is identical to the Bi equation in Ref. 
Again the solution must be the same: 


7 PyiX - aiPi, (A14) 

(where i was an atomic index). 


22 


P i(t + At) = g (At) (p i(t) - xAt'ypy^t)). (A15) 

Here g(At) = 1/y/l — 2c\At + C 2 AP 1 and cj and c 2 are expressions involving sums of 
products of molecular momenta (see Ref. 32). Next we consider the velocity relative 
to the center of mass. Dividing Eqs. f)A12[) and IA14I by rn ux and M t , respectively, 
and then subtracting the second from the first gives the equation for the relative 
velocity: 


d 

dt 


(Via - Vi) 


= 0, 


(A16) 


i.e. the relative velocity is not affected by the flow or the thermostat—this is of 
course the point of using molecular SLLOD equations and a molecular thermostat: 
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they translate each molecules rigidly, without rotation or deformation. The update- 
expression for an atomic velocity is therefore that which gives the same numerical 
difference as the update-expression for the center of mass velocity. 


Vi„(t + At) = Vi a (t) - Vi(f) + Vi (t + At) 

= v to (i) + (g(A t) - l)vi(i) - g(At)Atgv yl x. 


(All) 

(A1S) 


For the B2 equations of motion, the procedure is essentially the same: the velocity 
of the center of mass obeys the B 2 equation in Ref. 32 for the atomic case, with 
solution 


p i(t + At) = 


e — h/e 


h ( , . „ . 1 + h — e — h/e 


(1 - b)p 


= p i(t) + MjAvj (A19) 


where expressions for e(At), /3 and h are given in Ref. 33 (note that the sign of the 
second term in parentheses was incorrect in Ref. [321, their Eq. (22)). The notation 
Avj for the change in center of mass velocity of molecule i is introduced for use 
below. The relative velocity obeys: 


d_ 

dt 


(Via - Vi) = F ia /mi a 


Fj/ilTj n ; ; 0: a i; 


(A20) 


where & ia and a* are the (constant) accelerations. To update the atomic velocities 
we then have 


v ia (t + At) = v ia (t) + A Vi + At(a ia - a*), (A21) 

where An* is determined from Eq. (1A19(1 . 

A final technical point: while RUMD uses single precision floating point arith¬ 
metic in general, double precision is used for velocities in SLLOD algorithm. In this 
work, however, to keep the drift of the kinetic energy acceptably low for the small 
time step and strain rates used, double precision was used for the entire simulation 
code. 


41 











REFERENCES 


4 J. F. Branthaver, J. C. Pedersen, R. E. Robertson, J. J. Duvall, S. S. Kim, P. M. 
Harnsberger, T. Mill, E. K. Ensley, F. A. Barbour, and J. F. Schabron, Tech. Rep. 
SHRP-A-368, Strategic Highway Research Program (1993). 

2 D. A. Anderson, D. W. Christensen, H. U. Bahia, R. Dongre, M. G. Sharma, C. 
E. Antle, and J. Button, Tech. Rep. SHRP-A-369, Strategic Highway Research 
Program (1994). 

3 J. C. Petersen, Transportation Research Circular E-C140, Transportation Re¬ 
search Board (2009). 

4 S. E. Johnson, W. Y. Svrcek, and A. K. Mehrotra, Ind. Eng. Chem. Res. 26, 2290 
(1987). 

5 L.-I. Palade, P. Attane, and S. Carnaro, Rheol. Acta 39, 180 (2000). 

6 D. Sybilski, Mater. Struct. 26, 15 (1993). 

7 0. Ukwuoma and B. Ademodi, Fuel Processing Technology 60, 95 (1999). 

8 M. Garcia-Morales, P. Partal, F. J. Navarro, F. Martinez-Boza, C. Gallegos, N. 
Gonzalez, O. Gonzalez, and M. E. Munoz, Fuel 83, 31 (2004). 

9 P. Michalica, I. B. Kazatchkov, J. Stastna, and L. Zanzotto, Fuel 87, 3247 (2008). 
10 J. M. Dealy and J. Wang, Melt Rheology and its Applications in the Plastics 
Industry (Springer, New York, 2013). 
n O. C. Mullins, Annu. Rev. Anal. Chem. 4, 393 (2011). 

12 J. Murgich, J. Rodriguez M., and Y. Aray, Energy & Fuels 10, 68 (1996). 

13 J. H. Pacheco-Sanchez, I. P. Zaragoza, and J. M. Martinez-Magadan, Energy & 
Fuels 17, 1346 (2003). 

14 F. Frigerio and D. Molinari, Comp. Theor. Chem. 975, 76 (2011). 

15 P. Ungerer, D. Rigby, B. Leblanc, and M. Yiannourakou, Mol. Simulat. 40, 115 
(2014). 

16 L. Zhang and M. L. Greenheld, J. Chem. Phys. 127, 194502 (2007). 

17 J. S Hansen, C. A. Lemarchand, E. Nielsen, J. C. Dyre, and T. Schrpder, J. Chem. 
Phys. 138, 094508 (2013). 

18 C. A. Lemarchand, T. B. Schrpder, J. C. Dyre, and J. S. Hansen, J. Chem. Phys. 
139, 124506 (2013). 


42 



19 D. D. Li and M. L. Greenfield, Fuel 115, 347 (2014). 

20 P. Coussot and C. Ancey, Phys. Rev. E 59, 4445 (1999). 

21 R. G. Winkler, D. A. Fedosov, and G. Gompper, Curr.Opin. Colloid Interface Sci. 
19, 594 (2014). 

22 A. Vaccaro and G. Marrucci, J. Non-Newtonian Fluid Mech. 92, 261 (2000). 

23 D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids 
(Academic, London, 1990). 

24 B. D. Todd and P. J. Daivis, Molecular Simulation 33, 189 (2007). 

25 T. C. Le, B. D. Todd, P. J. Daivis, and A. Uhlherr, J. Chem. Phys. 131, 044902 
(2009). 

26 Z. Li, H. Djohari, and E. E. Dormidontova, J. Chem. Phys. 133, 184904 (2010). 

27 N.-T. Van-Oanh, C. Houriez, and B. Rousseau, Phys. Chem. Chem. Phys. 12, 930 
(2009). 

28 C02 emission reduction by exploitation of rolling resistance modelling of pave¬ 
ments, http://www.cooee-co2.dk/. 

29 C. A. Lemarchand, T. B. Schrpder, J. C. Dyre, and J. S. Hansen, J. Chem. Phys. 
141, 144308 (2014). 

30 A. Wiehe and K.S. Liang, Fluid Phase Equil. 117, 201 (1996). 

31 ASTM. 1995, Annual Book of Standards (American Society for Testing and Ma¬ 
terials, Philadelphia, 1995); method D-2007. 

32 G. Pan, J. F. Ely, C. McCabe, and D. J. Isbister, J. Chem. Phys. 122, 094114 
(2005). 

33 N. Bailey, L. Bphling, J. S. Hansen, T. Ingebrigsten, H. Larsen, C. A. Lemarchand, 
U. R. Pedersen, T. Schrpder, and A. A. Veldhorst, http://rumd.org. 

34 J. H. Irving and J. G. Kirkwood, J. Chem. Phys. 18, 817 (1950). 

35 R. B. Bird, C. F. Curtiss, R. C. Armstrong, and O. Hassager, Dynamics of Poly¬ 
meric Liquids (Wiley, New York, 1987). 

36 P. J. Farrington, C. J. Hawker, J. M. Frechet, and M. E. MacKay, Macromolecules 
31, 5043 (1998). 

37 C.-Y. Liu, J. He, R. Keunings, and C. Bailly, Macromolecules 39, 8867 (2006). 

38 X. Lu and U. Isacsson, Constr. Build. Mater. 14, 79 (2000). 


43 



39 J. G. Guan, M. Kariznovi, H. Nourozieh, and J. Abedi, J. Chem. Eng. Data 58, 
611 (2013). 

40 P. Partal, F. Martmez-Boza, B. Conde, and C. Gallegos, Fuel 78, 1 (1999). 

41 M. Masoori and M. L. Greenfield, J. Chem. Phys. 141, 124504 (2014). 

42 S. Bair, C. McCabe, and P. T. Cummings, Phys. Rev. Lett. 88, 058302 (2002). 
43 M. Van Gurp and J. Palmen, Rheol. Bull. 67, 5 (1998). 

44 J. Ge, B. D. Todd, G. Wu, and R. J. Sadus, Phys. Rev E 67, 061201 (2003). 

45 L. Separdar, N. P. Bailey, T. B. Schrpder, S. Davatolhagh, and J. C. Dyre, J. 
Chem. Phys. 138, 154505 (2013). 

46 J. T. Bosko, B. D. Todd, and R. J. Sadus, J. Chem. Phys. 121, 12050 (2004). 
47 M. Doi and S. F. Edwards, The theory of polymer dynamics (Clarendon, Oxford, 
1986). 

48 S. Forster, M. Konrad, and P. Lindner, Phys. Rev. Lett. 94, 017803 (2005). 

49 B. Aguilera-Mercado, C. Herdes, J. Murgich, and E. A. Muller, Energy & Fuels 
20, 327 (2006). 

50 J. L. S. Wales, The Application of Flow Birefringence to Rheological Studies of 
Polymer Melts (Delft LIniversity Press, Delft, 1976). 

51 R. Khare, J. de Pablo, and A. Yethiraj, J. Chem. Phys. 107, 6956 (1997). 

52 D. B. Genovese, Adv. Colloid Interface Sci. 171-172, 1 (2012). 

53 M. L. Mansfield and J. F. Douglas, J. Chem. Phys. 139, 044901 (2013). 

54 J. T. Bosko, B. D. Todd, and R. J. Sadus, J. Chem. Phys. 123, 034905 (2005). 
55 T. Indei, T. Koga, and F. Tanaka, Macromol. Rapid Commun. 26, 701 (2005). 
56 A. Tripathi, K. C. Tam, and G. H. McKinley, Macromolecules 39, 1981 (2006). 
57 T. Koga and F. Tanaka, Macromolecules 43, 3052 (2010). 

58 M. Kroger and R. Makhlouh, Phys. Rev. E 53, 2531 (1996). 

59 A. Milchev, J. Wittmer, and D. Landau, Eur. Phys. J. B 12, 241 (1999). 

60 D. Bedrov, G. D. Smith, and J. F. Douglas, Europhys. Lett. 59, 384 (2002). 

61 J. Padding, E. Boek, and W. Briels, J. Chem. Phys. 129, 074903 (2008). 

62 C. C. Huang, H. Xu, and J. Ryckaert, Europhys. Lett. 131, 58002 (2008). 

63 F. Zhang, D. J. Searles, D. J. Evans, J. S. den Toorn Hansen, and D. J. Isbister, 
J. Chem. Phys. Ill, 18 (1999). 


44 



